function PICCSTest(Phase,movements,noiseLevel)
%the size of the pic's resolution 
n =64;
N = n*n;

%set projection num
L = 180;

%measure matrix
matrix = getSampMatrix(n,L);
A = double(matrix);

%scan angle sqc for fbp
theta = 1:(180/L):180;

%dyna sample
cycleNum = 5;
cycleLen = floor(L/cycleNum);
halfCycleLen = floor(cycleLen/2);

% 
if(movements == 0 )
   MaxDis = 4;
else
   MaxDis = 2;
end

speed = MaxDis/halfCycleLen;
ObjSize = 4;

if halfCycleLen == 0
    disp('need more projection')
    return;
end

for i = 1:L
    j = mod(i,cycleLen);
    pos = j - halfCycleLen;
    dis = pos * speed;
    dis = floor(dis);
    a(j+1) = dis;
    if(movements == 0)
        samplePIC = SampPIC(30+dis,30+ObjSize+dis,30,30+ObjSize,n);
    else
        samplePIC = SampPIC(30-dis,30+ObjSize+dis,30-dis,30+ObjSize+dis,n);
    end

    Rtemp = radon(samplePIC,theta(i));
    R(:,i) = Rtemp;
    Ttheta(j+1,ceil(i/cycleLen)) = i;
    
    btemp = A * samplePIC(:);
   %-----------------add noise----
   noise = noiseLevel * randn(size(btemp));
   if(noise ~= 0)    
        btemp = btemp + noiseLevel * noise/norm(noise);
   end
	%-------------------------------
        proj(j+1,(ceil(i/cycleLen)-1)*n+1:(ceil(i/cycleLen))*n) = (i-1)*n+1:i*n;
        bcycle(j+1,(ceil(i/cycleLen)-1)*n+1:(ceil(i/cycleLen))*n) = btemp((i-1)*n+1:i*n);
    %b((i-1)*n+1:i*n,:)= btemp((i-1)*n+1:i*n,:);       
end

xpDyna = iradon(R,theta,'spline','Shepp-Logan',n);




%----------reconstruction-------------------------------------
%set lamda
lamda = 0;

%Phase = [30:32];


detVal =   bcycle(Phase,:)';
detVal = detVal(:);
PrjSeq = proj(Phase,:)';
PrjSeq = PrjSeq(:);
thetaSeq = Ttheta(Phase,:)';
thetaSeq = thetaSeq(:);

%CS reconstruction
cvx_begin
    variable xcs(N)
    %minimize( norm(xrec,1) )  
    minimize(lamda*norm(xcs - xpDyna(:),1)+ (1-lamda)*norm(xcs,1))
    %minimize (lamda*sum_square(xrec-xpDyna(:)) + (1-lamda)*norm(xrec,1))
    subject to
     %   bcycle(5,0*n+1:5*n)'  == A(proj(5,0*n+1:5*n),:) * xrec ;   
     detVal  == A(PrjSeq,:) * xcs ;   
     if(noise ~= 0)    
         norm(detVal - A(PrjSeq,:) * xcs)<noiseLevel ;  
     else
         detVal == A(PrjSeq,:) * xcs;
     end
cvx_end


%PICCS reconstruction
lamda = 0.2;
cvx_begin
    variable xpiccs(N)
    minimize(lamda*norm(xpiccs - xpDyna(:),1)+ (1-lamda)*norm(xpiccs,1))
    subject to
     if(noise ~= 0)    
         norm(detVal - A(PrjSeq,:) * xpiccs)<noiseLevel ;  
     else
         detVal == A(PrjSeq,:) * xpiccs;
     end  
cvx_end

%FBP reconstruction
bp= iradon(R(:,thetaSeq),theta(thetaSeq),64);

figure(1)
subplot(2,2,1)
showPIC(xpDyna(:),n,n);
title('Prior Image')
subplot(2,2,2);
showPIC(bp(:),n,n);
title('FBP reconstruction')
subplot(2,2,3)
showPIC(xcs,n,n);
title('CS reconstruction');
subplot (2,2,4)
showPIC(xpiccs,n,n);
title('PICCS reconstruction');


